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ABSTRACT 


Single-channel "pilot" manual control output in closed- 
tracking tasks is modeled in terms of linear discrete transfer 
functions which are parsimonious and guaranteed stable . The 
transfer functions are found by applying a modified super- 
position time series generation technique . A Levinson-Durbin 
algorithm is used to determine the filter which prewhi'tens the 
input and a projective ( least squares ) fit of pulse response 
estimates is used to guarantee indentif ied model stability. 
Results from two case studies are compared to previous 
findings , where the source of data are relatively short data 
records , approximately 25 seconds long . Time delay effects and 
pilot seasonalities are discussed and analyzed . It is 
concluded that sing le-channel time series controller modeling 
is feasible on short records , and that it is important for the 
analyst to determine a criterion for "best time domain fit" 
which allows association cf model parameter values , such as 
pure time delay, with actual physical and physiological 
constraints . The "purpose" of the modeling is thus paramount . 
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NOMENCLATURE 

numerator discrete polynomial in z 
coefficient of z _Jc in a(z ) 
denominator discrete polynomial in z 
coefficient of z“ k in b( z ) 
discrete pilot model pulse response 
error displayed to pilot at instant t 
coefficient of z~^ in g ( z ) 

pilot transfer function as a ratio of polynomials 
independent, identically distributed 
lag implying "kA" seconds 

pilot gain expressed in degrees per degree, 
total points available 

pilot input uncorrelated with y ( t ) in degrees at 
instant t 

white noise sequence (i.i.d. ) at instant t 
controlled element output signal in deg rees pitch 
angle at instant t 
sample interval ( seconds ) 

pilot output in degrees of elevator deflection at 
instant t 

number of sample times in pure time delay 

transformed frequency 

frequency 

prewhitening filter in z 
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1 . INTRODUCTION 


The key question of how the human being will be inser- 
ted in the control loop of complex processes remains an issue 
throughout our society (Rosenbrock, 1983), but nowhere is it 
more urgent than in flight control systems design and analysis 
(Harper , 1983). The fact that a pilot of a modern ai rcraf t is 
becoming" a sophisticated systems monitor (Rouse , 1983) in no 
way implies his demise as a controller (Rouse, 1980; Sheridan, 
1974), and a fundamental assumption in this work is that the 
interaction between man and machine should be understood much 
better than it is today (Palmer, 1983). 

Although describing function (McRuer , 1965 ) and optimal 
control (Kleinman, 1969-1 974 ) pilot models have been inge- 
niously used to provide insight into piloting strategy 
( Schmidt, 1979; Bacon, 1983; Hess , 1977), they are now supple- 
mented with pilot models derived from the emerging field of 
time series analysis. Time series modeling of pilot behavior 
offers tremendous potential for discerning key system charac- 
teristics and relationships , such as the actual effect of 
instabilities (Goto, 1974 ) , pilot stress (Shinners, 1974), or 
task effects (Agarwal, 1980 ) . 

The key questions in time series models involve not only 
the parsimony of parameters, well established by Breddermann et 
al (1978), but of identified model stability and the' model's 
practical application in analysis ( Baron, 1980). Shinner 
(1974) seriously discussed the closed-loop identification 
problem, but the manipulation of transfer functions in his 
fitting procedure contains no guarantee of final model stabi- 
lity. The primary purpose of this work is to present a theore- 
tically sound and relatively simple closed-loop fitting 
procedure, still based firmly in the common sense methods of 
Box and Jenkins (1976), which guarantees model stability 
without sacrificing model accuracy. 


2. MODEL 

The linear discrete closed-loop model structure is shown 
in Figure 1 . Each block represents a discrete pulse response 
sequence which, when convolved with the discrete input 
sequence, yields the discrete output sequence . Stable pulse 
sequences, even though inf inite in durations, eventually must 
decay for a stable system. When the pulse sequence is ex- 
pressed as a ratio of polynomials, stability is guaranteed if 
the denominator roots are less in magnitude than one. The goal 
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is to identify the pulse response sequence g p (z) and approxi- 
mate its discrete (z domain) transfer function from actual data 
sets |6(t)}, {y(t )} , and { W ( t ) } which are equispaced in time 
with their means removed. 


The assumptions are model linearity, time invariance, 
causality, uncorrelated inputs W(t) and R(t), and prewhitenable 
input W ( t ) ; that is, W(t) is a linear function of previous 
values plus a white noise "shock." Previous values are mathe- 
matically linked by the backward shift operator 


3. MODIFIED SUPERPOSITION TECHNIQUE 

First, every signal in Figure 1 is decomposed conceptually 
into a part linearly correlated with command disturbance W(t), 
the remainder uncorrelated with W(t). For example, output 
y(t) is the sum of Y^t), which is correlated with W(t), and of 
Y R (t), considered the effect of an additional unknown input 
R(t), termed "Remnant," uncorrelated with W(t). The pulse 
response to be found relates, for constant sampling interval 
"A" seconds, the linearly correlated pilot output 6 L (t) to the 
correlated error signal eL ( t ) ; that is, 

e L ( z )g p ( z ) = 6 l (z ) ( 1 ) 

This pulse response may be expressed as an infinite sequence or as 
a ratio of polynomials: 


Gp(z) = Kz T a(z)/b(z) = z 1 ( £ g k z ^ ) 

k = 0 


( 2 ) 
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I f the integer "k" in Equation (2) is allowed all values 
(_» < <+«>), then Equation (2) defines the discrete transfer 

function relating the z-transf orm of input sequence e L ( t ) to 
the z-transform of output sequence 6 L (t ) (Franklin and 
Powell , 1980, p. 1 5 ) . 

Although the signals ^L(t) and eL ( t ) are not directly 
available, they must be "generated" if loop closure effects are 
properly taken into account. To do this apply superposition to 


signals Y ( t ) and W(t) of Figure 1 : 

Y ( t ) = G 1 (z)W(t) + G 2 ( z ) R ( t ) (5) 

Y L ( t ) = G 1 (z)W(t) (6) 

where 

G 2 (z) = -G a (z)/[ 1-G a (z)G p (z) ] (7) 

G i ( z ) = Gp ( z )G 2 ( z ) ( 8 ) 


Since W(t) is prewhitenable (defined above) and uncorrelat- 
ed with R(t) the cross correlation identification technique of 
Box and Jenkins (details in Appendix) may be applied to find 
an estimate of the initial portion of the pulse response sequence 
gi(z), between y(t) and W(t). Then a-j(z) and b^(z) may be 
determined as shown in the section on model stability, such 
that 


G^ (z) = a-) (z.)/b} (z) (9 ) 

The essence of modified superposition is now to generate the time 
series Y L (t) using the autoregressive relation 

b 1 (z)Y L (t) = a 1 (z)W(t) (10) 

Where a^(z) and bi(z) are numerator and denominator polyno- 
mials, 

respectively, with the structure of Equations (3) and (4). The 
linearly correlated signal e^(t) is then generated from 

e L ( t ) =W ( t ) - Y L ( t ) (11) 

The above process is then repeated by reapplying super- 
position to obtain the following relation between Mt) and W(t): 


6(t) = G 3 ( z ) W ( t ) + G 4 (z)R(t) (12) 

6 L (t) = G 3 ( z ) W ( t ) (13) 
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The cross correlation identification (Appendix) applied to the 
sequence 6(t) and W(t) yields the initial segment of pulse 
response sequence 93(2)/ and the polynomials a3(z) and b3(z) may 
be determined (see next section) such that 


G3 ( z ) = a3(z)/b3(z) ( 14 ) 

Pilot output linearly correlated with W(t) is generated from the 
autoregressive relation 

b 3 ( z ) 6 L (t) = a3 ( z )W(t ) ( 15 ) 

Finally , the cross correlation technique (Appendix) is 
applied to 6 L (t) and e L (t) to find the initial segment of gp(k), 
defined by the' coefficient set {gp k r 0 5 k < N} , of the pilot model 
pulse response. Numerator and denominator polynomials are then 
found (see next section) which yields 

Gp( z ) = 6 l (z )/e L (z ) ( 16 ) 

No multiplication or divisions of transfer functions occurs 
throughtout the above procedure. 


4 . MODEL STABILITY 


As mentioned above, the pulse response sequence identifi- 
ed [ g 1 ( 2 ) , g3(z) and gp(z)] will be truncated at some finite 
lag "k H final task is to find a parsimonious numerator polyno- 
mial stable denominator polynomial which together are equivalent 
mathematically to the identified pulse response. These polyno- 
mials are chosen to have the structure shown in Equation ( 2 ), 
which is re-arranged into the following form: 


s . kmax , Ji, 

( 1+1 biZ _1 )( l g k z ) = K (1 +1 a k z K ) ( 17 ) 

i= 1 k=0 k= 1 

kmax >s -SL 
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Since the pulse response g k i s known for 0 <k <N, by equat- 
ing coefficients for the operator "z" at each exponential 
power "H" , o relationships may be found between numerator and deno- 
miator coefficients a^ and b^. Moreover, by equating coefficients 
for the operator z above power "s H , for which the right side of 
equation (15) vanishes, one obtains for every j > 0 

9s+ j + 9 s +j-1 b 1 + ••• + 9 s +j-s b s = 0 (18) 


The above relation exists for a finite but large number of 
H J >0", so projection theory (least squares) may be used to coef- 
ficients bk (0 <k -s). Bringing term M g( s +j)" to the other 
side of equation (18) and divided by "g(s + j) B one may write 

T T 

A[bi , b2, • . . , b s ) = [-1, ..., -1] (19) 

and the "j"th row 0 f A i s given by 

| g s+ j-1 , 9s+ j-2 , . . . , 9s+ j-s | j > 0 ( 20 ) 

g s+ j 9 3+ j g s+ j , by s 


The solution from linear algebra is 

T t — 1 T T 

[ fa 1' b 2 • ***' b s ^ = -(A A) A [ 1 , 1 ] (21) 

To provide a parsimonious denominator, the solution of 
Equation (21) is accepted for the lowest order "s" which has 
both a stable characteristic equation (i.e. roots less than 1.0 
in magnitude) and which yields a model pulse response similar 
in shape to the truncated pulse response identified from the 
data. Once a stable denominator is found the numerator a (z) 
and the gain K may be determined by once again matching coef- 
ficients in Equation (17); 

K = g Q (22) 


a k 


- K Uk 


i=k 

+ b i9k-i^ 


0 < k i i 


(23) 
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By defining error residual to be the actual output time 
series minus the pilot model output series at each sample 
instant, the gain K may be adjusted by a suitable minimization 
technique to minimize the error residual variance . 
Alternatively, it may be adjusted to provide a steady state 
response of unity when the input to the transfer function is a 
unity pulse train, a constraint recommended by Agarwal (1980). 

If a time delay "x " is to be included, the final form 
of Gp( z ) will be as shown in Equation (2), and the indices for 
the pulse responses in Equations ( 1 7 ) - ( 23 ) should be incre- 
mented by the integer "x " during identification (for example 
the gain K from Equation (2) equals g T identified from the 
data) . 


Validation tests may also be applied to the model . There 
are two types of tests : acceptability and statistical signifi- 
cance . Acceptability tests are common sense checks which com- 

pare model output series verses actual autocorrelation 
estimates from the data, autocorrelation of residuals for whi- 
teness properties , and checks for negligible cross -cor relation 
between the noise inputs . 

Statistical significance tests may be performed after 
acceptability tests indicate the model is reasonable . Chi- 
squared statistics are available from the w ( t ) and v ( t ) prewhi- 
tened series (discussed in the Appendix and shown in Figure 
13). Assuming one can safely neglect correlations beyond a lag 
of 20, for example, the statistics to be computed are, for 
“whiteness" of v(t) 

20 N 

(N-p ) l { l v(t-k) v( t) ( 24 ) 

k=1 (N-k ) t*k 

and, for uncorrelated w(t ) and v(t) 

20 . N 

(N-P) 1 { — l w ( t~k ) v ( t ) } (25) 

k s 1 (N-k) t=k 

p » order of a v (z ) filter 

N s total points in data set 

which should pass the chi-squared significance test for degrees 
of freedom ( 20-p ) and ( 20-1 -s-1 ) respectively (Box and Jenkins , 
1976, p.394). Failure of either significance test is evidence 
of a faulty assumption or a modeling inadequacy. 
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To summarize the modified superposition technique 

a) Find a finite pulse sequence relating y(t) and W ( t ) using 
cross correlation identification (Appendix ) . 

b) Determine a parsimonious , stable transfer function G •) ( I ) 
which is mathematically equivalent, in the least squares 
sense , to the sequence identified from the data g i ( z ) 

[ Equation ( 9 ) J . 

c ) Generate time series realizations { yj^ ( t ) } , (e^t)} us ing 
Equations (10) and (11). 

d) Find a finite pulse sequence, g3 (z ) , relating 6 ( t ) and W(t ) 
using cross correlation identification, and determine a 
stable transfer function 63(1) for this pulse response 
( Equation (14)). 

e ) Generate time realization 6 L ( t ) using Equation (15). 

f ) Find a finite sequence of the pulse response g p ( z ) , from 
6 L ( t ) and 6 L ( t ) using cross correlation identification, 
and fit a stable pilot model transfer function G p ( x ) to 
this pulse response (Equation (16)). 

g) Adjust K if desired and validate the model . 


5. PILOTED LABORATORY SIMULATION 

S ingle- channel "piloted" simulations in the Flight 

Simulation Laboratory at Purdue University were accomplished 
with a pilot performing pursuit tracking tasks us ing a single 
and double integrator (K/s and K/s^ respectively) controlled 
element dynamics . The task involved a command disturbance 

input of a random appearing forcing function, and a standard 
pursuit (McRuer , 1974) display using a CRT Monitor . Data sets 

were obtained at a 20 hertz sample rate and 500 points were 

used for modeling , providing a record length of only 25 seconds 
(although the data run itself exceeded 60 seconds ) . 

For the single- integrator controlled element many low- 
order transfer functions provided excellent "fits, " and the 
lowest order model is shown in Table 1. A "direct iden- 

tification" neglecting the closed-loop structure was also per- 
formed by merely fitting signals {6 ( t )} and { e ( t ) } , and a 
comparison of those results in Table 1 shows little variation 
in parameter values between direct and indirect identification 
in ' this case . This implies a small value for pilot injected 
noise relative to stick output (See Figure 1), a reasonable 
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deduction for a "simple" controlled element such as K/s . An a 
priori selected time delay of 0. 2 seconds yielded the lowest 
error residual variance and is consistent with previous results 
(Bredderman, 1976). 

A frequency response of the identified transfer func- 
tion is shown in Figures 2 and 3 where it is clear that a delay 
in series with a pure gain effectively describes pilot beha- 
vior. This is consistent with classical pilot modeling results 
(McRuer, 1974 ) . Since a conventional Bode interpretation and 
analysis using these frequency responses is not valid over all 
frequencies in discrete systems z-domain analysis , a transfor- 

om z to w ' was accomplished using 
10, p.114) 


(27) 

show the transformed frequency (uj 
response in the w' domain, where a conventional Bode interpre- 
tation is allowed. By comparing Figures 4 and 5 with Figures 2 
and 3, one can find no discernable difference between the 
responses over the frequency range of interest ( 0 to 25 rps ) . 

The time histories are shown in Figure 6. Only the first 
500 points (25 seconds ) were used to develop the model , -and the 
model output remains reasonably accurate beyond this time . This 
verifies stationar ity and avoids an overfit (Kashyap, 1976), 
which would be evidenced by increased error residual when the 
model is applied to data independent of model derivation ( in 
this case beyond 25 seconds ) . 

For the double- integrator controlled element a more 
complex transfer function was identified and is shown in Table 
2 for two values of a priori selected time delay (0.05 seconds 
and 0. 2 seconds ) . 

From the frequency response plot in Figure 7 there is some 
resonance near 2.0 Hz. The phase plots are shown in Figures 8 
and 9 for two different values of time delay (0.2 and 0.05 
seconds respectively) . The transformation to w ' domain yields 
no discernable difference from these responses and they are not 
shown. 


mation of variables f r 
(Franklin and Powell, 1 9 1 


w 1 


2 
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(z-1 
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tan 
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In contrast to control of a "simple" K/s , the "best fit 
(minimizing residual error variance) was obtained when time 
delay was set to 0.05 control for of K/s^ . The phase contribu- 
tion (from the poles and zeros) of the discrete transfer func- 
tion is apparent as time delay changes between 0.2 seconds and 
0.05 seconds, as may be seen by the phase plots of Figures 10 
and 11 in which the pure time delay has been removed from the 
discrete, transfer function. Selecting the larger pure time 
delay for the model exposes the considerable lead generation 
from the transfer function poles and zeros. This lead genera- 
tion is not as apparent when pure time delay is reduced for the 
"best time domain" fit, but the resulting 0.05 seconds might be 
judged too fast to associate with a lumped physiological delay 
for a human operator. A possible explanation is unmodeled 
pilot anticipation; that is, a possible anticipatory loop clo- 
sure not accounted for in Figure 1. 

Further evidence of this is provided in the time history 
for the best fitting model in Figure 12. Note that a seasonal 
pilot residual (where pilot output "leads" model output) occurs 
during some of the longer intervals of large slope. This could 
be caused by momentary anticipatory behavior arising from the 
"pursuit" display including commanded input, a factor not 
accounted for in a time invariant model. Thus in determining 
the "best" model using time series analysis, the purpose of the 
model must be given as much consideration as tests for "best 
fit." 


In summary for the K/s^ controlled element, an a priori 
time delay in series with a rate sensitive gain describes "pilot" 
behavior over his usable bandwidth, in agreement with classical 
results (McRuer, 1974 ). When pure time delay is not set a 
priori but allowed to vary in obtaining the "best time domain 
fit," the minimization of an error variance criterion results 
in a math model where the time delay is perhaps too small to be 
associated with physiological operator delays. This case is 
associated with a pursuit task in which the command as well as 
the plant output is displayed. 
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6. CONCLUSIONS 


A modified superposition technique was described for 
obtaining a parsimonious and stable discrete transfer function , 
along with statistical tests for model validation. Results 
provide evidence that the time series technique appears 
feasible to implement on "short" data records . The analyst 
needs, hQwever , to determine the criterion for a "best time 
domain fit" which allows association of parameter values , such 
as pure time delay, with actual physical and physiological 
constraints. Seasonalities in pilot residual , possibly caused 
by anticipatory behavior, we re observed as first noted by 
Shinners (1974), and are not well modeled with a time invariate 
model. 

Future work should concentrate on the full potential of 
these time series models for analyses, especially their ability 
to provide stable and accurate power spectral densities , and on 
their application to multi-channed closed-loop pilot modeling . 
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8. APPENDIX; Cross Correlation Identification 

(Box and Jenkins , 1976) 

Given the situation in Figure 13, the goal is to find the- 
pulse response relating Y ( t ) and W ( t ) , which is prewhi tenable 
by w w ( z ) . The prewhitening is accomplished by applying the 
Levinson-Durbin algorithm as given by Kay and Marple ( 1981 , pp . 
1388-1 389 ). By reversing the order of the blocks in the for- 


ward path of Figure 13 , and multiplying each signal at the 
summer by 

il -1 ( z ) , the following equation results ; 

w 

G ( z ) w ( t ) + $l w _1 ( z ) V ( t ) = 8 ( t ) (28) 

8 ( t ) = a w (z)y(t) (29) 

Now multiply Equation ( 28 ) by w ( t-k ) and take the expecta- 
tion, recalling that w ( t ) is uncorrelated by assumption 

with v( t ) : 

G ( z ) E [ w ( t )w ( t-k ) ] = E[ 8 (t )w(t-k) ] ( 30 ) 

By expanding G ( z ) using shift properties of z one obtains 
( 9o +< 9 1 2f 1 +92 z~ 2 + . • . ) E[w(t)w(t-k) ]= E[8 (t )w( t-k ) ] (31) 


Since w(t) is an independent, identically distributed sequence of 
random numbers with variance o w 2 t one obtains for every lag k 

9k°w 2 = E [ 8 ( t ) w ( t-k ) ] (32) 

k > 0 

Conventional estimation relations may now be used to estimate 
the terms in Equation ( 32 ) and solve for gk ; for example , from Box 
and Jenkins ( 1976, pp. 32-33 ) one obtains 


, N N 

g k t - 1 w(t)w(t)| = [— I 8 ( t ) w ( t-k ) } (33 ) 

N t=1 N-k t=k 

which detemines the pulse response sequence estimate g k- 
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Table 1 Discrete Transfer Function Identification 
Results for Controlled Element K/s (K * 1) 


K "^l+a.z" 1 ) 

Model Structure G (z) = ; — 

P Otbz' 1 ) 

Signal to noise ratio =50 K 

N ■ 500 points a » 0.05 seconds t = 4 (0.2 seconds) 


Parameter 

H 

■7TIIT ■ 

Direct 

Identification 

K P 

0.64 

0.69 

K** 

P 

0.79 

0.72 

a l 

0.71 

0.69 

b ! 

0.32 

0.37 


* Gain which minimizes error residual variance 
** Gain yields steady state step response of unity 
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Table 2 Modified Superposition Identification Results 
for Controlled Element K/s2 (K = 1) 


Kz^l+a.z' 1 ) 

Model Structure G (z) = 3 3 3— 

p (l+b 1 z~ 1 +b 2 z"^+b ; j2”' 3 ) 

Signal to noise ratio = 30 

N « 500 a = 0.05 seconds t = 4 (0.2 seconds) 


Parameter 

t = 0.05 sec 

i = 0.2 seconds 

K 

k p 

0.03 

0.89 

V* 

0.033 

1.22 

a i 

10.9 

-0.67 

“i 

- 1.42 

-1.41 

b 2 

0.91 

0.88 

“3 

- 0.1 

-0.06 

Roots 

0.14 

0.64 ±j 0.57 

0.08 

0.67 ±j 0.57 


* Gain which minimizes error residual variance 
** Gain yields steady state step response of unity 
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Figure 4 W' Respond Magnitude: Controlled Element K/i 
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Figure 5 W' Response Phase: Controlled Element Kh 
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Figure 6 Model Output vs. Pilot Output: Controlled Element 6C/S 
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DEGREES 


Figure 8 Manual Controller Frequency Response Phase: 
K/s 2 , r * 0.2 seconds 
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